program principal

implicit none

    interface
        subroutine bounding_surfaces_read_in(BOUNDING, n, m, xmin, xmax, ymin, ymax, nonvalue)
            integer :: n, m
            real :: xmin, xmax, ymin, ymax, nonvalue
            real, allocatable :: BOUNDING(:,:,:)
        end subroutine
    end interface

    interface
        subroutine readzmapfilesMUD(Numsurf,MUDDRAPE,SPACE, MUDDRAPEREAL, filepath)
            real, allocatable :: MUDDRAPE(:,:)
            real, allocatable :: MUDDRAPEREAL(:,:)
            integer :: Numsurf,SPACE
            character :: filepath*1000
        end subroutine
    end interface
    
    interface
        subroutine writegridfaciesfileMULTZ(GRIDFACIESMULTZ,nblayers,n,m, filepath6)
            integer :: m,n,nblayers
            character :: filepath6*1000 
            real, allocatable :: GRIDFACIESMULTZ(:,:,:)
        end subroutine
    end interface
    
    interface
        subroutine writegridfaciesfileMULTY(GRIDFACIESMULTY,nblayers,n,m, filepath6)
            integer :: m,n,nblayers
            character :: filepath6*1000
            real, allocatable :: GRIDFACIESMULTY(:,:,:)
        end subroutine
    end interface
    
    interface
        subroutine writegridfaciesfileMULTX(GRIDFACIESMULTX,nblayers,n,m, filepath6)
            integer :: m,n,nblayers
            character :: filepath6*1000
            real, allocatable :: GRIDFACIESMULTX(:,:,:)
        end subroutine
    end interface
    
    interface
        subroutine tenpercentZ(GRIDFACIESMULTZ)
            real, allocatable :: GRIDFACIESMULTZ(:,:,:)
        end subroutine
    end interface
    
    interface
        subroutine tenpercentY(GRIDFACIESMULTY, dimi, dimj)
            integer :: dimi, dimj
            real, allocatable :: GRIDFACIESMULTY(:,:,:)
        end subroutine
    end interface
    
    interface
        subroutine tenpercentX(GRIDFACIESMULTX, dimi, dimj)
            integer :: dimi, dimj
            real, allocatable :: GRIDFACIESMULTX(:,:,:)
        end subroutine
    end interface
    
   interface 
    subroutine write_zmap_file(BOUNDING,temp_name,n,m,k,xmin,xmax,ymin,ymax)
    character :: temp_name*1000 
    integer :: n,m,k
    real :: xmin, xmax, ymin, ymax
    real, allocatable :: BOUNDING(:,:,:)
    end subroutine
   end interface
    
character :: temp_name*1000, filepath*1000, filepath2*1000, filepath3*1000, filepath4*1000, filepath5*1000, filepath6*1000    
integer:: Nbsurf,multi,multj, a, b, c, d, g 
real, allocatable :: MUDDRAPE(:,:)
real, allocatable :: MUDDRAPEREAL(:,:)
real, allocatable :: BOUNDING(:,:,:)
real, allocatable :: GRIDFACIESMULTZ(:,:,:)
real, allocatable :: GRIDFACIESMULTY(:,:,:)
real, allocatable :: MUDDRAPETEMP(:,:,:)
real, allocatable :: GRIDFACIESMULTX(:,:,:)
integer:: n,m,i,nbzones,j,k,nblayers,numlayer,nbcontrol,Numsurf,dimi,dimj
real:: Nonvalue,limits_surface(6),resolxgrid,resolygrid, xmin, xmax, ymin, ymax

call bounding_surfaces_read_in(BOUNDING, n, m, xmin, xmax, ymin, ymax, nonvalue)

!open(unit=8, file="mud coverage.txt")  
!close(unit=8, STATUS='DELETE')

write(*,*)"Enter the file path (including file name) for the Heterogeneity" 
write(*,*)"Surface information"
    read(*,*) filepath5
  
write(*,*)"Enter the file path (no file name necessary) for Heterogeneity Surfaces"
    read(*,*) filepath
    
write(*,*)"Enter the file path (including file name)for Grid information file"
    read(*,*) filepath4
    
write(*,*)"Enter the file path (including filename) Top Grid Layers Information file"
    read(*,*) filepath2
    
write(*,*)"Enter the file path (including filename) Base Grid Layers Information"
    read(*,*) filepath3
    
write(*,*)"Enter the output destination path for the grid property files"   
    read(*,*) filepath6 
     
open(unit=1,file=filepath5)
    read(1,*) Nonvalue
    read(1,*) limits_surface(1:6)
    read(1,*) n,m
    read(1,*) Nbsurf
close(unit=1)

open(unit=2, file=filepath4)
    read(2,*) nbzones
    read(2,*) nblayers
    read(2,*) dimi,dimj

multi=(n-1)/dimi
multj=(m-1)/dimj
!----------------------------------------------------------------------------------------------------------------------------------------------------
allocate(GRIDFACIESMULTZ(nblayers,dimj,dimi))
do k=1,nblayers
    do i=1,dimi
    do j=1,dimj
        GRIDFACIESMULTZ(k,j,i)=1
    end do
    end do
end do
  
print*, "Calculating MULTZ..............."

do k=1,nbzones

read(2,*) numlayer

    if (allocated(MUDDRAPE)) then;
    deallocate(MUDDRAPE)
    end if
    if (allocated(MUDDRAPEREAL)) then;
    deallocate(MUDDRAPEREAL)
    end if
    allocate(MUDDRAPE(n,m))
    allocate(MUDDRAPEREAL(n,m))

    !!!!!!!!!!!!!!!!!
    call readzmapfilesMUD(k,MUDDRAPE,n,MUDDRAPEREAL, filepath)
    
    do i=1,dimi
    do j=1,dimj
    a=(multi*i-1)
    b=(multj*j-1) 
    c=(multi*i)
    d=(multj*j) 
    if(a<1.or.b<1.or.c<1.or.d<1)then;
    GRIDFACIESMULTZ(numlayer,dimj-j+1,i)=1
    else;
        nbcontrol=MUDDRAPE(a,b)*100+MUDDRAPE(a,d)*1000+MUDDRAPE(c,b)*1+MUDDRAPE(c,d)*10
    end if
        if(nbcontrol==1001.or.nbcontrol==110.or.nbcontrol==1.or.nbcontrol==10.or.nbcontrol==100.or.nbcontrol==1000.or.nbcontrol==0) then;
        GRIDFACIESMULTZ(numlayer,dimj-j+1,i)=1
        else;
        GRIDFACIESMULTZ(numlayer,dimj-j+1,i)=0
        end if
    end do
    end do 

end do

close(unit=2)

!call tenpercentZ(GRIDFACIESMULTZ)

call writegridfaciesfileMULTZ(GRIDFACIESMULTZ,nblayers,dimj,dimi, filepath6)
deallocate(GRIDFACIESMULTZ)
!---------------------------------------------------------------------------------------------------------------------------------------------------------------
allocate(GRIDFACIESMULTY(nblayers,dimj,dimi))

do k=1,nblayers
    do i=1,dimi
    do j=1,dimj
        GRIDFACIESMULTY(k,j,i)=1
    end do
    end do
end do

print*, "Calculating MULTY..............."
  
open(unit=6, file=filepath2)

do g=1,1
    read(6,*) numlayer
  do k=1,nbzones

    if (allocated(MUDDRAPE)) then;
    deallocate(MUDDRAPE)
    end if
    if (allocated(MUDDRAPEREAL)) then;
    deallocate(MUDDRAPEREAL)
    end if
    allocate(MUDDRAPE(n,m))
    allocate(MUDDRAPEREAL(n,m))        

    call readzmapfilesMUD(k,MUDDRAPE,n,MUDDRAPEREAL, filepath)

    do i=1,dimi
    do j=1,dimj
    if(MUDDRAPEREAL(i,j)==BOUNDING(2,i,j))then;
    a=(multi*i-1)
    b=(multj*j-1) 
    c=(multi*i)
    d=(multj*j) 
    if(a<1.or.b<1.or.c<1.or.d<1)then;
    GRIDFACIESMULTY(numlayer,dimj-j+1,i)=GRIDFACIESMULTY(numlayer,dimj-j+1,i)
    else;
        nbcontrol=MUDDRAPE(a,b)*100+MUDDRAPE(a,d)*1000+MUDDRAPE(c,b)*1+MUDDRAPE(c,d)*10
    end if
        if(nbcontrol==1001.or.nbcontrol==110.or.nbcontrol==1.or.nbcontrol==10.or.nbcontrol==100.or.nbcontrol==1000.or.nbcontrol==0) then;
        GRIDFACIESMULTY(numlayer,dimj-j+1,i)=GRIDFACIESMULTY(numlayer,dimj-j+1,i)
        else;
        GRIDFACIESMULTY(numlayer,dimj-j+1,i)=0
        end if
    end if
    end do
   end do
    
  end do
end do

close(unit=6)

open(unit=6, file=filepath3)

do g=1,1
    read(6,*) numlayer
  do k=1,nbzones

    if (allocated(MUDDRAPE)) then;
    deallocate(MUDDRAPE)
    end if
    if (allocated(MUDDRAPEREAL)) then;
    deallocate(MUDDRAPEREAL)
    end if
    allocate(MUDDRAPE(n,m))
    allocate(MUDDRAPEREAL(n,m))        
    
    call readzmapfilesMUD(k,MUDDRAPE,n,MUDDRAPEREAL, filepath)
    
    do i=1,dimi
    do j=1,dimj
    if(MUDDRAPEREAL(i,j)>=BOUNDING(1,i,j))then;
    a=(multi*i-1)
    b=(multj*j-1) 
    c=(multi*i)
    d=(multj*j) 
    if(a<1.or.b<1.or.c<1.or.d<1)then;
    GRIDFACIESMULTY(numlayer,dimj-j+1,i)=GRIDFACIESMULTY(numlayer,dimj-j+1,i)
    else;
        nbcontrol=MUDDRAPE(a,b)*100+MUDDRAPE(a,d)*1000+MUDDRAPE(c,b)*1+MUDDRAPE(c,d)*10
    end if
        if(nbcontrol==1001.or.nbcontrol==110.or.nbcontrol==1.or.nbcontrol==10.or.nbcontrol==100.or.nbcontrol==1000.or.nbcontrol==0) then;
        GRIDFACIESMULTY(numlayer,dimj-j+1,i)=GRIDFACIESMULTY(numlayer,dimj-j+1,i) 
        else;
        GRIDFACIESMULTY(numlayer,dimj-j+1,i)=0
        end if
    end if
    end do
   end do
    
  end do
end do

close(unit=6)

!call tenpercentY(GRIDFACIESMULTY, dimi, dimj)

call writegridfaciesfileMULTY(GRIDFACIESMULTY,nblayers,dimj,dimi, filepath6)
deallocate(GRIDFACIESMULTY)
!-------------------------------------------------------------------------------------------------------------------------------------------------
allocate(GRIDFACIESMULTX(nblayers,dimj,dimi))

do k=1,nblayers
    do i=1,dimi
    do j=1,dimj
        GRIDFACIESMULTX(k,j,i)=1
    end do
    end do
end do

print*, "Calculating MULTX..............."

open(unit=6, file=filepath2)

do g=1,1
    read(6,*) numlayer
  do k=1,nbzones

    if (allocated(MUDDRAPE)) then;
    deallocate(MUDDRAPE)
    end if
    if (allocated(MUDDRAPEREAL)) then;
    deallocate(MUDDRAPEREAL)
    end if
    allocate(MUDDRAPE(n,m))
    allocate(MUDDRAPEREAL(n,m))        

    call readzmapfilesMUD(k,MUDDRAPE,n,MUDDRAPEREAL, filepath)

    do i=1,dimi
    do j=1,dimj
    if(MUDDRAPEREAL(i,j)==BOUNDING(2,i,j))then;
    a=(multi*i-1)
    b=(multj*j-1) 
    c=(multi*i)
    d=(multj*j) 
    if(a<1.or.b<1.or.c<1.or.d<1)then;
    GRIDFACIESMULTX(numlayer,dimj-j+1,i)=GRIDFACIESMULTX(numlayer,dimj-j+1,i)
    else;
        nbcontrol=MUDDRAPE(a,b)*100+MUDDRAPE(a,d)*1000+MUDDRAPE(c,b)*1+MUDDRAPE(c,d)*10
    end if
        if(nbcontrol==1001.or.nbcontrol==110.or.nbcontrol==1.or.nbcontrol==10.or.nbcontrol==100.or.nbcontrol==1000.or.nbcontrol==0) then;
        GRIDFACIESMULTX(numlayer,dimj-j+1,i)=GRIDFACIESMULTX(numlayer,dimj-j+1,i)
        else;
        GRIDFACIESMULTX(numlayer,dimj-j+1,i)=0
        end if
    end if
    end do
   end do
    
  end do
end do

close(unit=6)

open(unit=6, file=filepath3)

do g=1,1
    read(6,*) numlayer
  do k=1,nbzones

    if (allocated(MUDDRAPE)) then;
    deallocate(MUDDRAPE)
    end if
    if (allocated(MUDDRAPEREAL)) then;
    deallocate(MUDDRAPEREAL)
    end if
    allocate(MUDDRAPE(n,m))
    allocate(MUDDRAPEREAL(n,m))        
    
    call readzmapfilesMUD(k,MUDDRAPE,n,MUDDRAPEREAL, filepath)
    
    do i=1,dimi
    do j=1,dimj
    if(MUDDRAPEREAL(i,j)>=BOUNDING(1,i,j))then;
    a=(multi*i-1)
    b=(multj*j-1) 
    c=(multi*i)
    d=(multj*j) 
    if(a<1.or.b<1.or.c<1.or.d<1)then;
    GRIDFACIESMULTX(numlayer,dimj-j+1,i)=GRIDFACIESMULTX(numlayer,dimj-j+1,i)
    else;
        nbcontrol=MUDDRAPE(a,b)*100+MUDDRAPE(a,d)*1000+MUDDRAPE(c,b)*1+MUDDRAPE(c,d)*10
    end if
        if(nbcontrol==1001.or.nbcontrol==110.or.nbcontrol==1.or.nbcontrol==10.or.nbcontrol==100.or.nbcontrol==1000.or.nbcontrol==0) then;
        GRIDFACIESMULTX(numlayer,dimj-j+1,i)=GRIDFACIESMULTX(numlayer,dimj-j+1,i) 
        else;
        GRIDFACIESMULTX(numlayer,dimj-j+1,i)=0
        end if
    end if
    end do
   end do
    
  end do
end do

close(unit=6)

!call tenpercentX(GRIDFACIESMULTX, dimi, dimj)

call writegridfaciesfileMULTX(GRIDFACIESMULTX,nblayers,dimj,dimi, filepath6)
deallocate(GRIDFACIESMULTX)
!-------------------------------------------------------------------------------------------------------------------------------------------------
close(unit=3, STATUS='DELETE')  
print*, "Press any key to exit"
read(*,*)

stop
end program
!-------------------------------------------------------------------------------------------------------------------------------------------------
subroutine readzmapfilesMUD(Numsurf,MUDDRAPE,SPACE,MUDDRAPEREAL, filepath)

integer:: n,m,Numsurf,ibis,nfile,mfile,i,j,k,SPACE
real, allocatable:: MUDDRAPE(:,:)
real, allocatable:: MUDDRAPEREAL(:,:)
integer:: dummy1,lines,rest,h,raw,mode, numsurfCUBEmax
character:: name*50,ie*5,id*5, filepath*1000
character:: dummy2*10
real :: numb,xminfile,xmaxfile,yminfile,ymaxfile,Nonvaluefile
real :: passagefive(5),passagefour(4),passagethree(3),passagetwo(2)
real :: passageone
real:: passage(SPACE,5)

call listing(1000)

    open(unit=3,file="listing.txt")
        do ibis=1,Numsurf
            read(3,*) ie
        end do
    close(unit=3)

name="Heterogeneity_Surface"//trim(ie)//".zmap+"

open(unit=4,file=trim(filepath)//"Heterogeneity_Surface"//trim(ie)//".zmap+")
    read(4,*) dummy2
    read(4,*) dummy2
    read(4,*) dummy2
    read(4,*) dummy2
    read(4,*) dummy2
    read(4,*) dummy2
    read(4,*) numb, Nonvaluefile
    read(4,*) mfile, nfile, xminfile, xmaxfile, yminfile, ymaxfile
    read(4,*) dummy2
    read(4,*) dummy2

    lines=mfile/5+1
    rest=mfile-(lines-1)*5
    if (rest==0) then;
        lines=lines-1
    end if
    
    do i=1,nfile
        h=1
        do j=1,(lines-1)
        read(4,*) passagefive
        passage(h,1:5)=passagefive
        h=h+1
        end do

        if (rest==1) then;
          read (4,*) passageone
          passage(h,1)=passageone
        else if (rest==2) then;
          read (4,*) passagetwo
          passage(h,1:2)=passagetwo
        else if (rest==3) then;
          read (4,*) passagethree
          passage(h,1:3)=passagethree
        else if (rest==4) then;
          read (4,*) passagefour
          passage(h,1:4)=passagefour
        else;
          read (4,*) passagefive
          passage(h,1:5)=passagefive
        end if
        
        do k=1,mfile
          raw=k/5
          mode=k-raw*5
          if (mode==0) then;
            mode=5
            raw=raw-1
          end if
          MUDDRAPE(i,k)=passage(raw+1,mode)          
        end do
     end do
 close(unit=4)
 
  do i=1,nfile
 do j=1,mfile
     MUDDRAPEREAL(i,j)=MUDDRAPE(i,j)
  end do
 end do
 
 do i=1,nfile
 do j=1,mfile
     if (MUDDRAPE(i,j)-1<Nonvaluefile) then;
     MUDDRAPE(i,j)=0
     else;
     MUDDRAPE(i,j)=1
     end if
 end do
 end do
 
 end subroutine readzmapfilesMUD
!---------------------------------------------------------------------------------------------------------------------------------------------------
subroutine writegridfaciesfileMULTZ(GRIDFACIESMULTZ,nblayers,n,m, filepath6)

character :: name*50, filepath6*1000
integer :: m,n,i,j,nblayers 
real, allocatable :: GRIDFACIESMULTZ(:,:,:)

    open(unit=5,file=trim(filepath6)//"MULTZ.GRDECL")
    write(5,"(A5)") "TRANZ"
    do k=1,nblayers
    do i=1,n
    do j=1,m
           write(5,*) GRIDFACIESMULTZ(k,i,j)
    end do
    end do
    end do
    write(5,"(A1)") "/"
    close(unit=5)
   
end subroutine writegridfaciesfileMULTZ
!---------------------------------------------------------------------------------------------------------------------------------------------------
subroutine writegridfaciesfileMULTY(GRIDFACIESMULTY,nblayers,n,m, filepath6)

character :: name*50, filepath6*1000
integer :: m,n,i,j,nblayers 
real, allocatable :: GRIDFACIESMULTY(:,:,:)

    open(unit=8,file=trim(filepath6)//"MULTY.GRDECL")
    write(8,"(A5)") "TRANY"
    do k=1,nblayers
    do i=1,n
    do j=1,m
           write(8,*) GRIDFACIESMULTY(k,i,j)
    end do
    end do
    end do
    write(8,"(A1)") "/"
    close(unit=8)
   
end subroutine writegridfaciesfileMULTY
!-----------------------------------------------------------------------------------------------------------------------------------------------------
subroutine writegridfaciesfileMULTX(GRIDFACIESMULTX,nblayers,n,m, filepath6)

character :: name*50, filepath6*1000
integer :: m,n,i,j,nblayers 
real, allocatable :: GRIDFACIESMULTX(:,:,:)

    open(unit=9,file=trim(filepath6)//"MULTX.GRDECL")
    write(9,"(A5)") "TRANX"
    do k=1,nblayers
    do i=1,n
    do j=1,m
           write(9,*) GRIDFACIESMULTX(k,i,j)
    end do
    end do
    end do
    write(9,"(A1)") "/"
    close(unit=9)
   
end subroutine writegridfaciesfileMULTX
!-----------------------------------------------------------------------------------------------------------------------------------------------------
subroutine listing(SPACE)

 integer :: SPACE, i

 open(unit=3, file="listing.txt")
  do i=1,SPACE
   write(3,*) i
  end do
 close(unit=3)

 return
end subroutine listing
!-------------------------------------------------------------------------------------------------------------------------------
subroutine bounding_surfaces_read_in(BOUNDING, n, m, xmin, xmax, ymin, ymax, nonvalue)

  interface
    subroutine array_allocation (temp_name,TEMP,n,m, nonvalue,xmin, xmax, ymin, ymax, passage)
    character :: temp_name*1000 
    integer :: n,m
    real :: xmin, xmax, ymin, ymax, nonvalue
    real, allocatable :: TEMP(:,:,:)
    real, allocatable :: passage(:,:)
    end subroutine
   end interface
    
    character :: temp_name*1000
    integer :: k, n, m, SPACE
    parameter (SPACE=6100)
    real :: nonvalue, xmin, xmax, ymin, ymax   
    real,allocatable :: BOUNDING(:,:,:)
    real, allocatable :: TEMP(:,:,:)
    real, allocatable :: passage(:,:)
    
    allocate (passage(SPACE,5))
    allocate (TEMP(2,SPACE,SPACE))
    
    write(*,*)"Enter the file path (including file name) of the base bounding surface:"
    read(*,*) temp_name
    call array_allocation (temp_name,TEMP,n,m, nonvalue,xmin, xmax, ymin, ymax, passage)  
    allocate (BOUNDING(2,(m*2),(m*2)))
    k=1
        do i=1,n
            do j=1,m
               BOUNDING(k,i,j)=TEMP(1,i,j)
            end do
        end do

    write(*,*)"Enter the file path (including file name) of the top bounding surface:"
    read(*,*) temp_name
    call array_allocation (temp_name,TEMP,n,m, nonvalue,xmin, xmax, ymin, ymax, passage)
    k=2
        do i=1,n
            do j=1,m
               BOUNDING(k,i,j)=TEMP(1,i,j) 
            end do
        end do
     
     deallocate (passage)
     deallocate (TEMP)    
     
return
end subroutine bounding_surfaces_read_in
!-------------------------------------------------------------------------------------------------------------------------------
subroutine array_allocation (temp_name,TEMP,n,m, nonvalue,xmin, xmax, ymin, ymax, passage)

    character :: temp_name*1000
    character :: dummy1*10, dummy2*10, dummy3*10, dummy4*10
    character :: dummy5*10, dummy6*10, dummy7*10, dummy8*10
    integer :: m, n, lines, rest, i, j, h, raw, mode, k, SPACE
    parameter (SPACE=6100)
    real :: numb, nonvalue, xmin, xmax, ymin, ymax, depth_min, depth_max
    real :: passagefive(5), passagefour(4), passagethree(3), passagetwo(2)
    real :: passageone
    real, allocatable :: passage(:,:)
    real, allocatable :: TEMP(:,:,:)  
    
   open(unit=1, file=temp_name)
	 read(1,*) dummy1
	 read(1,*) dummy2
	 read(1,*) dummy3
	 read(1,*) dummy4
	 read(1,*) dummy5
	 read(1,*) dummy6
	 read(1,*) numb, nonvalue
	 read(1,*) m, n, xmin, xmax, ymin, ymax
	 read(1,*) dummy7
     read(1,*) dummy8
     
     do i=1,SPACE
       do j=1,SPACE
          TEMP(1,i,j)=nonvalue
          TEMP(2,i,j)=nonvalue
       end do
    end do
      
     lines=m/5+1
     rest=m-(lines-1)*5
	if (rest==0) then;
    	lines=lines-1
    end if
    
	 do i=1,n
    	h=1
        do j=1,(lines-1)
           read(1,*) passagefive
           passage(h,1:5)=passagefive
           h=h+1
    	end do

		if (rest==1) then;
          read (1,*) passageone
          passage(h,1)=passageone
        else if (rest==2) then;
          read (1,*) passagetwo
          passage(h,1:2)=passagetwo
        else if (rest==3) then;
          read (1,*) passagethree
          passage(h,1:3)=passagethree
        else if (rest==4) then;
          read (1,*) passagefour
          passage(h,1:4)=passagefour
        else;
          read (1,*) passagefive
          passage(h,1:5)=passagefive
		end if
        
		do k=1,m
          raw=k/5
          mode=k-raw*5
          if (mode==0) then;
            mode=5
            raw=raw-1
          end if
          TEMP(1,i,k)=passage(raw+1,mode)          
         end do
     end do
    close(unit=1) 
              
return
end subroutine array_allocation
!-------------------------------------------------------------------------------------------------------------------------------
subroutine write_zmap_file(TEMP,temp_name,n,m,k,xmin,xmax,ymin,ymax)

    character :: temp_name*28
    character :: path*97
    integer :: n, m, i, k
    real :: xmin, xmax, ymin, ymax
    real, allocatable :: TEMP(:,:,:)
          
    path="C:\Users\gg809\Desktop\"

     open(unit=1,file=path//temp_name)
        write(1,*)'!'
        write(1,*)'!     ZIMS FILE NAME :  Clinoform_Surface.txt'
        write(1,*)'!     FORMATTED FILE CREATION DATE: xx/xx/xxxx'
        write(1,*)'!     FORMATTED FILE CREATION TIME: xx:xx'
        write(1,*)'!'
        write(1,1)'@Grid HEADER, GRID, 5'
1       format(A21)
        write(1,*)'	15,-9999.000,  ,	4 ,1.0000'
        write(1,"(I6,A,I5,A,1F12.4,A,1X,1F12.4,A,1F12.4,A,1F12.4)") m,',',n,',',xmin,',',xmax,',',ymin,',',ymax
        write(1,*)'0.0000     0.0000     0.0000'
	    write(1,2)'@'
2       format(A1)
        do i = 1, n 
            write(1,"(1F10.4,1X,1F10.4,1X,1F10.4,1X,1F10.4,1X,1F10.4)") TEMP(k,i,1:m)
         end do
    close(unit=1)    

return
end subroutine write_zmap_file
!-------------------------------------------------------------------------------------------------------------------------------------------------------
subroutine tenpercentZ(GRIDFACIESMULTZ)

integer :: dimj,dimi,i,j,nblayers, c, count, numlayer, nbzones,d
real, allocatable :: GRIDFACIESMULTZ(:,:,:)


open(unit=2, file="layers.txt")
    read(2,*) nbzones
    read(2,*) nblayers
    read(2,*) dimi,dimj

c=0
count=0    
do k=1,nbzones

read(2,*) numlayer

    do i=1,dimi
    do j=1,dimj
    if(GRIDFACIESMULTZ(numlayer,j,i)==0)then;
      c=c+1
     end if   
    end do
    end do

count=(c/2)
d=0
  do i=1,dimi
   do j=1,dimj
    if (GRIDFACIESMULTZ(numlayer,j,i)==0)then;
    GRIDFACIESMULTZ(numlayer,j,i)=0
    d=d+1
    end if
    if(d>count)then;
    GRIDFACIESMULTZ(numlayer,j,i)=1
    d=d+1
    end if
    end do
    end do
    
d=0
c=0
count=0       
end do
close (unit=2)

return
end subroutine tenpercentZ
!------------------------------------------------------------------------------------------------------------------------------------------------------------------
subroutine tenpercentY(GRIDFACIESMULTY, dimi, dimj)

integer :: dimj,dimi,i,j,nblayers, c, count, numlayer, nbzones
real, allocatable :: GRIDFACIESMULTY(:,:,:)

c=0
count=0  
open(unit=6, file="layersTOP.txt")

do g=1,2
    read(6,*) numlayer

    do i=1,dimi
    do j=1,dimj
    if(GRIDFACIESMULTY(numlayer,j,i)==0)then;
      c=c+1
     end if   
    end do
    end do

count=c
c=0
do while(c<(count-2))
  do i=1,dimi
   do j=1,dimj
    if (GRIDFACIESMULTY(numlayer,j,i)==0)then;
    c=c+1
    if(MOD((c/2),2)==0)then;
    GRIDFACIESMULTY(numlayer,j,i)=1
    end if
    end if
   end do
  end do      
end do

end do
close (unit=6)

return
end subroutine tenpercentY
!------------------------------------------------------------------------------------------------------------------------------------------------------------------
subroutine tenpercentX(GRIDFACIESMULTX, dimi, dimj)

integer :: dimj,dimi,i,j,nblayers, c, count, numlayer, nbzones
real, allocatable :: GRIDFACIESMULTX(:,:,:)

c=0
count=0  
open(unit=6, file="layersTOP.txt")

do g=1,2
    read(6,*) numlayer

    do i=1,dimi
    do j=1,dimj
    if(GRIDFACIESMULTX(numlayer,j,i)==0)then;
      c=c+1
     end if   
    end do
    end do

count=c
c=0
do while(c<(count-2))
  do i=1,dimi
   do j=1,dimj
    if (GRIDFACIESMULTX(numlayer,j,i)==0)then;
    c=c+1
    if(MOD((c/2),2)==0)then;
    GRIDFACIESMULTX(numlayer,j,i)=1
    end if
    end if
   end do
  end do      
end do

end do
close (unit=6)

return
end subroutine tenpercentX